function fp = perturbation(sc, auxdata)

Re = auxdata.Re;
mu = auxdata.mu;
J2 = auxdata.J2;

p = sc(1);
f = sc(2);
g = sc(3);
h = sc(4);
k = sc(5);
L = sc(6);

[~,e,i,Ome,ome,nu] = lowThrustMee2Coe(p,f,g,h,k,L);

[r,~] = orb2rv(p,e,i,Ome,ome,nu,L,nu+ome,Ome+ome,mu);
r = norm(r);
fr = -3*mu*J2*Re^2/(2*r^4) * ( 1 - 12*(h*sin(L) - k*cos(L))^2 / (1+h^2+k^2)^2 );
ft = -12*mu*J2*Re^2/(r^4) * ( (h*sin(L) - k*cos(L))*(h*cos(L)+k*sin(L)) / (1+h^2+k^2)^2 );
fh = -6*mu*J2*Re^2/(r^4) * ( (1-h^2-k^2)*(h*sin(L) - k*cos(L)) / (1+h^2+k^2)^2 );

fp = [ft;fr;fh];
